Movable Type Scripts

Distances and bearings between points on an ellipsoidal-model earth

For the benefit of the terminally obsessive (as well as the genuinely needy, of course), Thaddeus Vincenty devised formulæ for calculating geodesic distances between a pair of lati­tude/longi­tude points on the earth’s surface, using an accurate ellipsoidal model of the earth.

When I looked up the ref­erences (‘Direct and Inverse Solutions of Geodesics on the Ellipsoid with applica­tion of nested equa­tions’), I discovered to my surprise that while the mathematics is utterly beyond me, it is actually quite simple to program.

Vincenty’s solution for the distance between points on an ellipsoidal earth model is accurate to within 0.5 mm distance (!), 0.000015″ bearing, on the ellipsoid being used. Calcula­tions based on a spherical earth model, such as the (much simpler) Haversine, are accurate to around 0.3% – which is still good enough for many (most?) purposes, of course.

Vincenty’s inverse solution can fail on nearly antipodal points. Testing with Geographic­Lib test data, I’ve found this can happen with distances greater than 19,936 km, or within around 75 km of the antipodal point.

Vincenty wrote this at a time when computing resources were expen­sive, and made his solu­tions very computa­tionally efficient – the distance calcula­tion takes around 5 – 10 micro­seconds (hence around 100,000 – 200,000 per second) using Chrome on an aging Core i5 PC.

Distance/bearing between two points (inverse solution)

Vincenty’s paper presents a set of formulæ rather than a single program; this shows those formulæ pulled together as they are used in the script:

Where:

• s is the geodesic distance along the surface of the ellipsoid (in the same units as a & b)
• α₁ is the initial bearing, or forward azimuth
• α₂ is the final bearing (in direction p₁→p₂)

Note: Vincenty observes that eqn. (18) becomes indeterminate over equatorial lines (since cos²α → 0); in this case, set cos2σ_m to 0 and the result is computed correctly. He also points out that the formula may have no solu­tion between two nearly antipodal points; an itera­tion limit traps this case (Vincenty says “this will occur when λ, as computed by eqn. (11), is greater than π in absolute value”, but in a small number of cases, the formula will not converge even when λ stays below π).

Note: some implementations of Vincenty’s formula inefficiently use a large number of trig functions; Vincenty devised this solu­tion with an eye for efficiency in implementa­tion, and this one uses just one each of sin, cos, sqrt, and atan2 for each itera­tion – only 3 or 4 itera­tions are generally required. Pythagorean trigonometric identities sin²φ + cos²φ = 1 and 1 + tan²φ = sec²φ are used to reduce the number of trig functions invoked.

JavaScript

const L = λ2 - λ1; // L = difference in longitude, U = reduced latitude, defined by tan U = (1-f)·tanφ.
const tanU1 = (1-f) * Math.tan(φ1), cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1 * cosU1;
const tanU2 = (1-f) * Math.tan(φ2), cosU2 = 1 / Math.sqrt((1 + tanU2*tanU2)), sinU2 = tanU2 * cosU2;

let λ = L, sinλ = null, cosλ = null;    // λ = difference in longitude on an auxiliary sphere
let σ = null, sinσ = null, cosσ = null; // σ = angular distance P₁ P₂ on the sphere
let cos2σₘ = null;                      // σₘ = angular distance on the sphere from the equator to the midpoint of the line
let cosSqα = null;                      // α = azimuth of the geodesic at the equator

let λʹ = null;
do {
    sinλ = Math.sin(λ);
    cosλ = Math.cos(λ);
    const sinSqσ = (cosU2*sinλ) * (cosU2*sinλ) + (cosU1*sinU2-sinU1*cosU2*cosλ)**2);
    sinσ = Math.sqrt(sinSqσ);
    cosσ = sinU1*sinU2 + cosU1*cosU2*cosλ;
    σ = Math.atan2(sinσ, cosσ);
    const sinα = cosU1 * cosU2 * sinλ / sinσ;
    cosSqα = 1 - sinα*sinα;
    cos2σₘ = cosσ - 2*sinU1*sinU2/cosSqα;
    const C = f/16*cosSqα*(4+f*(4-3*cosSqα));
    λʹ = λ;
    λ = L + (1-C) * f * sinα * (σ + C*sinσ*(cos2σₘ+C*cosσ*(-1+2*cos2σₘ*cos2σₘ)));
} while (Math.abs(λ-λʹ) > 1e-12);

const uSq = cosSqα * (a*a - b*b) / (b*b);
const A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
const B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
const Δσ = B*sinσ*(cos2σₘ+B/4*(cosσ*(-1+2*cos2σₘ*cos2σₘ)-B/6*cos2σₘ*(-3+4*sinσ*sinσ)*(-3+4*cos2σₘ*cos2σₘ)));

const s = b*A*(σ-Δσ); // s = length of the geodesic

const α1 = Math.atan2(cosU2*sinλ,  cosU1*sinU2-sinU1*cosU2*cosλ); // initial bearing
const α2 = Math.atan2(cosU1*sinλ, -sinU1*cosU2+cosU1*sinU2*cosλ); // final bearing

See below for the complete code.

Destination given distance & bearing from start point (direct solution)

Again, this presents Vincenty’s formulæ arranged to be close to how they are used in the script.

Where:

• φ₂, λ₂ is destination point
• α₂ is final bearing (in direction p₁→p₂)

JavaScript

const sinα1 = Math.sin(α1);
const cosα1 = Math.cos(α1);

const tanU1 = (1-f) * Math.tan(φ1), cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1 * cosU1;
const σ1 = Math.atan2(tanU1, cosα1); // σ1 = angular distance on the sphere from the equator to P1
const sinα = cosU1 * sinα1;          // α = azimuth of the geodesic at the equator
const cosSqα = 1 - sinα*sinα;
const uSq = cosSqα * (a*a - b*b) / (b*b);
const A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
const B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));

let σ = s / (b*A), sinσ = null, cosσ = null; // σ = angular distance P₁ P₂ on the sphere
let cos2σₘ = null; // σₘ = angular distance on the sphere from the equator to the midpoint of the line

let σʹ = null;
do {
    cos2σₘ = Math.cos(2*σ1 + σ);
    sinσ = Math.sin(σ);
    cosσ = Math.cos(σ);
    const Δσ = B*sinσ*(cos2σₘ+B/4*(cosσ*(-1+2*cos2σₘ*cos2σₘ)-B/6*cos2σₘ*(-3+4*sinσ*sinσ)*(-3+4*cos2σₘ*cos2σₘ)));
    σʹ = σ;
    σ = s / (b*A) + Δσ;
} while (Math.abs(σ-σʹ) > 1e-12);

const x = sinU1*sinσ - cosU1*cosσ*cosα1;
const φ2 = Math.atan2(sinU1*cosσ + cosU1*sinσ*cosα1, (1-f)*Math.sqrt(sinα*sinα + x*x));
const λ = Math.atan2(sinσ*sinα1, cosU1*cosσ - sinU1*sinσ*cosα1);
const C = f/16*cosSqα*(4+f*(4-3*cosSqα));
const L = λ - (1-C) * f * sinα * (σ + C*sinσ*(cos2σₘ+C*cosσ*(-1+2*cos2σₘ*cos2σₘ)));
const λ2 = λ1 + L;

const α2 = Math.atan2(sinα, -x); // final bearing

See below for the complete code.

Datums

While is it wonderful to be able to calculate distances to millimetre accuracy, there are certain caveats which it is worth being aware of if you are looking for accuracy better than around ±1 metre.

Millimetre precision is stretching the limits of geodetic survey techniques: very few points will have been surveyed to such precision (it’s easy to provide lati­tude/longi­tude ref­erences to 6 decimals of degrees; less easy to make that precision meaningful rather than spurious).

Further, the familiar WGS-84 (World Geodetic System 1984) has no ‘physical realisa­tion’ – it is not tied to geodetic groundsta­tions, just to satellites – and is defined to be accurate to no better than around ±1 metre, irrespective of how good your GPS unit is [actually, the ‘G-series’ are a little more accurate, but there are no official transforma­tion param­eters]. A metre is good enough for most purposes, but a long way from Vincenty’s millimetre accuracy.

Then, when working to an accuracy of better than ±1 metre, plate tectonic movements become significant. Depending on the datum it’s defined in, any given point will most likely have shifted many millimetres from where it was last year.

Simplifying somewhat (well, actually, a lot!), in order to address requirements for greater precision, the ITRS was developed to give greater accuracy than WGS84, with ‘epoch’-dependant ITRF ‘realisa­tions’, where the lati­tude/longi­tude coordinate of a posi­tion will vary over time.

Various ‘static’ reference frames are also defined for various continents – NAD83 for North America, ETRS89 for Europe, GDA94 for Australia, etc – within which lati­tude/longi­tude coordinates remain essentially fixed (at least to centimetre or so accuracy, major earthquake events excepted). These reference frames have ‘epoch’-dependant mappings to ITRF datums. (Due to plate tectonic continental drift, ETRS89 shifts against ITRF by about 25mm/year; GDA94 by around 80mm/year).

And for any latitude/longitude coordinates dating from earlier datums (such as NAD27, ED50, NTF, AGD66, and any number of others), even metre precision is exceeding their limits.

Also, it must be remembered that these distances are calculated on the surface of the reference ellipsoid, which will differ (to varying degrees) from the earth geoid, and even more from the earth’s physical surface. The geoid is equivalent to the local mean-sea-level (or its equivalent in land-locked locales), and is determined by the local force of gravity; it varies in height from the ellipsoid by around ±100 metres. The earth’s surface varies much more; if you happen to be located in Colorado, 2 kilometres above msl, distances will be 0.03% greater than on the ellipsoid surface; +1 metre over just 30 kilometres.

So for millimetre accuracy to be relevant, you need to know a fair bit about the reference frames, datums, and epochs of your coordinates, and other param­eters of your measure­ments. Caveat emptor!

The ubiquitous WGS-84 is a geocentric datum, based on an ellipsoid with:

Latitude/longitude points can be converted between different datums by converting them to geocentric cartesian coordinates, applying a 7-param­eter Helmert transforma­tion, and converting back to a geodetic coordinate in the new datum.

Should even millimetre accuracy not be sufficient for you, Charles Karney has improved on Vincenty’s method with a method which has errors in nanometers (and always converges on antipodal points) – Geodesics on an ellipsoid of revolu­tion, 2011, Algorithms for geodesics, 2012; though Karney’s method is a lot less straight­forward (someday I may get around to an implementa­tion).

Live examples

... to verify the code is correct!

Enter the co-ordinates into the text boxes to try out the calculations. A variety of formats are accepted, principally:

• deg-min-sec suffixed with N/S/E/W (e.g. 40°44′55″N, 73 59 11W), or
• signed decimal degrees without compass direction, where negative indicates west/south (e.g. 40.7486, -73.9864):

Test case (from Geoscience Australia), using WGS-84: